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*> ■ Abstract 
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<^ • We present a realistic expanding source model with nine parameters that are 

necessary and sufficient to describe the main physics occuring during hydro- 
dynamical freezeout of the excited hadronic matter produced in relativistic 
heavy- ion collisions. As a first test of the model, we compare it to data 
from central Si + Au collisions at pi a t>A4 = 14.6 GeV/c measured in exper- 
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iment E-802 at the AGS. An overall x 2 P er degree of freedom of 1.055 is 
achieved for a fit to 1416 data points involving invariant ir + , ir~ , K + , and 
K~ one-particle multiplicity distributions and ir + and K + two-particle cor- 
relations. The 99%-confidence region of parameter space is identified, leading 
to one-dimensional error estimates on the nine fitted parameters and other 
calculated physical quantities. Three of the most important results are the 
freezeout temperature, longitudinal proper time, and baryon density along the 
symmetry axis. For these we find values of 92.9 ± 4.4 MeV, 8.2 ± 2.2 fm/c, 

and 0.0222 toiolt fm ~ 3 ' respectively 
PACS: 25.75.-q, 21.65. +f, 24.10. Jv, 24.10.Nz 
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I. INTRODUCTION 



It is a widely accepted theory that if nuclear matter attains a high enough energy density, 
it will undergo a phase transition from normal hadronic matter into a quark-gluon plasma 
(QGP) fl|-§|]. Since the discovery of such a QGP would represent a significant advancement 
in the fundamental understanding of nuclear interactions, there are a number of relativistic 
heavy-ion experiments both currently running and being planned which hope to test this 
theory. Unfortunately, if a QGP is formed in the laboratory, its quick expansion and cooling 
will cause it to transition back into normal hadronic matter long before anything can be 
detected. Thus, any signals for the prior existence of a QGP will necessarily be subtle and 
indirect. 

In order to work backwards from the final observed state of the detected hadrons to an 
earlier state which may or may not have included a QGP, it is necessary to use a reliable 
transport model. One approach which has been quite successful in the past is to treat the 
expanding nuclear matter as a hydrodynamical fluid. This fluid is very hot and dense imme- 
diately after the collision, but with time it expands and cools. When some criterion is met 
(e.g., falling below a certain temperature or density), it is assumed that the fluid "freezes 
out" and becomes a collection of non-interacting, free-streaming hadrons. The freezeout 
hypersurface is thus some three-dimensional surface which separates hydrodynamically in- 
teracting nuclear fluid from free-streaming hadrons. According to this picture, when these 
hadrons are observed in detectors, their distributions and correlations contain information 
about the temperature, expansion velocities, chemical potentials, size, and shape of the fluid 
during freezeout. 

The purpose of this paper is to present a physically reasonable parametrization of the 
freezeout process, and then to find the best values for the freezeout parameters by comparing 
theoretical distributions and correlations to experimental data through a minimization of x 2 . 
This approach is somewhat different from a standard nuclear hydrodynamical approach, in 
which some equation of state must be assumed in order to determine how the fluid evolves 
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from its initial condition to its final freezeout ||. One problem with standard nuclear 
hydrodynamics is that the formidable computations involved make a minimization of x 2 
impractical, so even when the agreement with experiment is quite good, one can never be 
sure that the best point in the infinite-dimensional space of all possible initial conditions, 
equations of state, and freezeout criteria has been found. Our more limited goal is to tackle 
just the problem of determining the properties of the system during freezeout. 

We begin by reviewing the Wigner-function formulation of hydrodynamical freezeout 
and defining nine parameters that are necessary and sufficient to properly describe the gross 
properties of the source during freezeout. Although in Sec. II. B we use the language of 
hydrodynamical evolution (e.g., rarefaction waves and cooling) to motivate our approach, it 
should be noted that our calculations are actually concerned only with freezeout — not with 
the hydrodynamical evolution which might have led to it. Section II. D then includes a short 
explanation of how resonance decays are taken into consideration. Once the model is defined, 
Sec. Ill outlines our general method for constructing x 2 , determining the goodness of the 
fit, and estimating uncertainties in the model parameters. With these tools in hand, we 
compare our nine-parameter model to data from central Si + Au collisions at pi a b/^4 = 14.6 
GeV/c, measured in experiment E-802 at the Brookhaven Alternating Gradient Synchrotron 
(AGS) p-g]. The 1416 data points used consist of invariant tt + , it", K + , and K~ one- 
particle multiplicity distribution measurements as well as tt + and K + two-particle correlation 
measurements. To our knowledge, this paper represents the first attempt to simultaneously 
find the best fit to one-particle distribution and three-dimensional two-particle correlation 
data with a single expanding source model. We found that the fits converged rapidly and 
consistently, yielding an overall x 2 P er degree of freedom of 1.055. 



II. DETAILS OF THE MODEL 
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A. Wigner Function Formulation 



The Wigner function for particles of type a with spin J a coming directly from a hydro- 
dynamical system involving a sharp three-dimensional freezeout hypersurface is |J 



^dir/ ^ _2J a + l p-n(x) 



a 1 ' P) (2tt) 3 exp{\p-u(x)-fi a (x)]/T(x)}Tl' 1 J 

where the — (+) sign is for bosons (fermions). The quantities u M (x), T(x), fi a (x), and 

n„(x) = [ d 3 a^x')5^(x-x') (2) 



s 

denote the local hydrodynamical flow velocity, temperature, chemical potential, and normal- 
pointing freezeout hypersurface element, respectively. Throughout the paper we use units 
in which h = c = k = 1, where ft is Planck's constant divided by 2ir, c is the speed of light, 
and k is the Boltzmann constant (except in the figures and tables, where we reinsert c). 
Integrating the direct Wigner function over spacetime generates the Cooper-Frye formula 



for the one-particle distribution [10 



^f(p) = / *xS? (x,p) = 2 -^±± [ d 3 a, P " , m • (3) 

J (27r) a Jt. exj)[{p-u{x) — fi a (x))/T{x)\ =f 1 

The subscript £ on the integral denotes the limits to the hypersurface for a finite-sized 
system. Because the observed particles are on mass shell, P^ ir depends only on the three- 
vector p rather than on the four-vector p. 

In addition to particles coming directly from the freezeout surface, there are also some 
which come from the decay of resonances. The total Wigner function for particle a is then 
comprised of two parts: 

S a (x,p) = S* 1T (x,p) + S ies -+ a (x,p) , (4) 

where the second term is determined by the direct Wigner functions of the contributing 
resonances (see Sec. II. D). The total observed multiplicity distribution for particle a is 

Pa(p) = / d 4 xS a (x,p) . (5) 



The correlation function for two particles of type a with momenta pi and P2 can similarly 
be expressed in terms of the Wigner function [0-II5] 



C a ( q ,K)-l±X a Pa(pi)Pa(p2) , (6) 

where the deviation from unity of the parameter X a measures the amount of coherent pro- 
duction of particles of type a. Although the on-shell momenta of the two particles is com- 
pletely specified by the six momentum components in K = |(pi + P2) arid q = Pi — P2, 
it is nevertheless notationally convenient to make the full off-shell four-vector definitions 
K = \(Vi + P2) and q = pi — p%- Since both the one-particle distribution and the two- 
particle correlation function are determined by the Wigner function, we need only find a 
suitable parametrization of this function in order to compare a hydrodynamical model to 
these data. 

B. Definition of the Model Parameters 

Our model is applicable to nearly central collisions of ultrarelativistic nuclei. For large 
sets of many nearly central collisions, the data should be azimuthally symmetric, so we 
assume azimuthal symmetry in our model. Immediately after the collision, we assume the 
formation of a hot, dense source which moves with some velocity v s = tanh?/ s relative to the 
lab while it expands and cools in its own rest frame. If the incoming nuclei are relativistic 
enough in the source frame, their strong Lorentz contraction makes their thickness in the 
beam (z) direction negligible, so it should be a good approximation to assume that the 
collision took place on a single plane at t — z — 0. Assuming also that the longitudinal flow 
velocities subsequently imparted to each bit of the nuclear fluid remain constant throughout 



the expansion, these velocities take the simple form first suggested in [|T^,|T^], namely 



P z (rj) = j = tanh?7 , (7) 

where r] = tanh _1 (2;/t) is the spacetime rapidity of that bit of fluid in the source frame. We 
will show shortly that this flow profile leads to a longitudinally boost-invariant local energy. 



Unlike in the longitudinal direction, there is no initial motion in the transverse direc- 
tion. After the collision, however, rarefaction waves work their way radially inward, causing 
the matter to accelerate transversally outward. The resulting three-dimensional expansion 
causes the fluid to cool until eventually a low enough temperature is reached so that the 
matter effectively stops interacting and "freezes out." We consider a model in which both 
the temperature T and chemical potential are constant at freezeout. For the latter, we define 

fl>a = B a fl h + S a fl B + I a p\ . (8) 

Here B a , S a , and I a are the baryon, strangeness, and isospin numbers of particle type a, 
while /it,, p s , and pi are the corresponding chemical potentials. 

Although there are many possible ways to parametrize the radial flow at freezeout, the 
actual profile chosen may not be nearly as important as the average transverse velocity of 
the profile |18[ . Recent hydrodynamical studies have obtained transverse flow profiles which 
are relatively linear in p = \lx 2 + y 2 out to a certain radius, outside of which they drop off 
quickly |T9|-f2T||. For simplicity, we assume a linear profile, and to preserve boost-invariance, 
we follow pi|-pB[] by defining it to be independent of z and t in the longitudinally comoving 



frame of the source. In other words, we parametrize the total flow velocity of the system in 
the source frame by 

u^{x) = 7p( cosh?/, /3 p cos0, /3 p sin0, sinh?/ J , (9) 



where 7 P = 1/y 1 — /3 P , with 

Mp) = * (£) • (io) 

Here R is the maximum transverse radius of the source and v t is the magnitude of the 
transverse velocity of the fluid at p = R. Note that f3 p is the flow velocity in the longitudinally 
comoving frame, but that the transverse component of the total flow velocity in the source 
frame is (3 P / cosh 77. 

That the flow profile is in fact boost-invariant can be most easily seen by first rewriting 
the source-frame particle four-momentum in the form 
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P IJr = {m t cosh(?/ p - y s ), p ± , m t sinh(?/p - y s ) ) , (11) 

where mt = \JE' 2 — p z ' 2 is the "transverse mass," y p = tanh -1 (p z /E) is the rapidity of the 
particle in the lab frame, and is the transverse momentum two-vector. Throughout 
this paper, we will use the subscript _L to denote the vector made from the two transverse 
components of a four-vector. Note that the quantity 

p-u(x) = j p m t cosh [(y p - y 8 ) - 77] - p_l-u_l(p) (12) 

depends on the rapidity of the particle and the spacetime rapidity of the source only through 
their difference. Since boosting to a frame moving with longitudinal velocity U relative to 
the source frame can be done by subtracting tanh _1 ([7) from both (y p — y s ) and r/, the 
difference of these quantities is boost invariant. 

In keeping with the boost-invariant profile, we assume that freezeout along the p = 
symmetry axis of the source occurs at a constant proper time [17]. Due to transverse expan- 



sion effects, however, freezeout may occur sooner for matter with p 7^ 0. These assumptions 
are incorporated into the following equation describing the freezeout hypersurface: 

t 2 -z 2 



Tf = const. , (13) 



l + a t (p/R) 2 

where a t parametrizes the radial behavior of the freezeout process. At a given constant slice 
in z, for — 1 < a t < 0, freezeout proceeds radially from outside to inside. For example, 



freezeout for the z = slice begins on the outside (p = R) at time t — t\ — t { ^1 + a t 
and continues until the inside (p = 0) freezes out at time t — ti — Tf. The case a t > 
corresponds to the less-likely possibility of freezeout proceeding radially from the inside to 
the outside, while a t = represents a freezeout which occurs at the same time for all points 
with a given z. The temporal duration of freezeout for the z = slice is just given by 



At(z = 0) = \t 2 - ti| = r f 1 - VT+at ■ (14) 

To derive the prefactor p-n(x) for the spacelike hypersurface defined by Eq. (0), it is 
most convenient to use the spacelike variables rj, x, and y. We have [^l[ 



j3 , , dX»dX a dXP 

d a^(x) = e^p— — — — — drjdxdy , (15) 



where the coordinate vector on the hypersurface is given by 



X» = ^T { ^l + a t {p/R) 2 cosh V , x, y, r f yj 1 + a t {p/R) 2 sinh^ J . (16) 

Doing the algebra, we find 

p-d 3 a(x) = r f |\/l + a t (p/R) 2 m t cosh [(y p - y s ) - rj] - a t p ± -x±Tf/ R 2 ^ drjdxdy , (17) 

where xj_ = (x, y). Again, since the only (y p — y s ) and rj dependencies come through the 
difference [(y p — y s ) — rj], the above expression is longitudinally boost invariant. 

Since we are interested in realistic finite systems, it is necessary to put some spacelike 
limits on the hypersurface. We have already mentioned the maximum transverse radius 
R; we also assume the existence of a maximum longitudinal radius z 3 = n sinh r] Q which is 
achieved by the source at time £3 = n cosh r]o- Because colliding nuclei have more matter in 
the center (p = 0) than on the outside, we take our source to be spheroidal in p and r\ (as 
opposed to cylindrical, for example). In other words, on the hypersurface E, the spacelike 
coordinates p and 77 satisfy the inequality 

£ + (18) 

R l T]q a 

Since in the above equation, rj appears alone and not in combination with (y p — y B ), these 



limits break boost invariance. Using Eq. (|2|) and the definition r = \Jt 2 — z 2 , we obtain in 
the source frame 

p-n(x) = < m t cosh [(y p - y B ) - rj\ '"' P± X± ' ' 



R 2 ^/l + a t (p/R) 2 

xs(t- r f ^i + a t (pAR) 2 ) e (1 - { P /Rf - (v/m) 2 ) . (19) 

The freezeout hypersurface as a function of t, p, and z shown in Fig. 1 corresponds to 
R = 8.0 fm, Tf = 8.2 fm/c, 770 = 1.47, and a t = —0.86, which are parameters that we 
determine in Sec. IV.B for the reaction considered here As mentioned before, since 



a t is negative, freezeout begins on the outside at p — R and works its way in, reaching 
the center at t — t 2 = r { . Figure 2 shows an illustration of the hydrodynamical fluid for 
seven different instances in source-frame time. The inner surfaces at each time are actually 
freezing out, while the outer end caps are the boundaries of fluid which will freeze out later. 
Since the end caps are not yet on the freezeout hypersurface, their exact shapes are not 
actually determined by our model and are shown only to illustrate the finite nature of the 
source. By time t = Tf, freezeout has worked its way to the center, so for later times the 
source becomes two separated receding fireballs of continually decreasing size. 

For symmetric projectile-target collisions, the rapidity y s of the source frame relative to 
the lab is given simply by the average of the projectile and target rapidities. For asymmetric 
collisions, however, the precise value of y s depends on how many "participant" nucleons in 
each nucleus collide to form the hydrodynamical source. The center of mass of the source is 
just the center of mass of the incoming participants rather than the total center of mass of 
all of the incoming nucleons (participants + spectators). Given the masses of the incoming 
nuclei, y s could be estimated either purely on geometrical grounds or it could be treated as 
another variable parameter to be fit to data. For the asymmetric Si + Au collisions studied 
in this paper, we choose the latter approach. 

Although at first we considered separate incoherence parameters for pions and for kaons, 
we subsequently found that very good fits could be obtained by setting X K — 1 and allowing 
only A,,- to vary as a parameter. Moreover, two of the chemical potentials, /i s and /!;, can be 
determined from the remaining nine parameters by imposing the constraints that the total 
strangeness of the sum of all particles in the source vanishes and that the total isospin per 
baryon is the same as that of the participants before the collision (see next subsection). The 
nine adjustable parameters of our model can be grouped in the following way: T, fih/T, 
and Ajr describe intrinsic properties of the fluid; R, v t , and a t describe transverse aspects 
of the freezeout; and y s , rjo, and Tf describe longitudinal and boost-invariant aspects of the 
freezeout. 
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C. Constraint Equations 



The total number of particles of type a which freeze out on the hypersurface is given by 



N? = J ^^(P) , (20) 

where P^ r is defined by Eq. @. Given the numbers of each particle species, the total 
strangeness and isospin per baryon of the system are simply 

Stot(/4, Mi) = S aN^{fJ, s , /ii) 
a 

Aot (Ats, Mi) = EcJaN^Mh) , 21 , 

where we have suppressed all parameter dependencies except those of /x s and /ii. The sum 
in a is over all mesons with masses below 900 MeV and all baryons with masses below 1410 
MeV. For particles of a given mass, all of the isospin, baryon, and strangeness states are 
considered separately since the different chemical potential of each (see Eq. (|8|)) leads to a 
different value of N^ ir . 

The initial isospin per baryon of the system depends on the number of participant protons 
and neutrons from each nucleus. We define the target proton number, nucleon number, and 
number of participants as Z t3I , A tax , and B tar , respectively. Making similar definitions for 
the projectile and noting that each proton (neutron) has isospin 1/2 (—1/2), we find the 
total isospin of the incoming participants: 

B B 

I\o — JTZ [-^proj — (^4proj ~~ -Zproj)] + [^tar ~~ (^tar — -^tar)] • (22) 

^^iproj ^^tar 



To get the isospin per baryon of the participants, we simply divide the above equation by 
B = -Bproj + B tax . The quantities B pro ^/B and B tar /B can be determined by equating the 
incoming target and projectile momenta in the participant center-of-mass frame. Explicitly, 

m Nj B proj sinh(y pro j - y s ) = m Nj B tar sinh(?/ s - y tax ) , (23) 

where y pTO j and y ta ,r are the initial projectile and target rapidities and is the nucleon 
mass. Using this relation, it is easy to show that the initial isospin per baryon of the system 
is given by 
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1 (ZprojMproj) sinh(?/ s - y tar ) + (Z tar /A tar ) sinh(y, 



proj 




(24) 



B o 



2 sinh(?/ s - y t3X ) + sinh(y proj - y a ) 



In general, the initial isospin per baryon depends upon the parameter y s , but for symmetric 
collisions (or any collision in which Z pTO j/A pTO j = Z tar /A tar ), I/B\ is independent of y s . 
We are now ready to write down explicitly the constraint equations for /x s and ix\. 



For a given set of the nine remaining parameters, these equations allow one to find unique 
values of /x s and /i; as well as their derivatives (e.g., dfi a /dT). 

Notice that the baryon number of this model reflects only the participant baryons and is 
not constrained to be the same as the total baryon number of the two incoming nuclei. The 
rationale behind this is that in many collisions there are "spectator" nucleons whose evolu- 
tion is not well described by a hydrodynamically expanding source. Since these spectators 
may nonetheless end up in detectors, this model works best when it is used to fit only data 
from produced particles such as mesons and antiprotons. Of course, once the parameters 
have been determined from a fit to mesons, for example, it is always possible to compare 
the proton distribution predicted by the model with the data to get an idea of how many of 
the measured protons are "participants" and how many are "spectators." 

It is significant that even if only meson data are considered, it is still possible to determine 
the baryon chemical potential /ib from the relative abundances of K +, s and K~'s. For any 
system which has a positive baryon number, there are more A's and S's produced than 
A's and S's. This leads to a net negative strangeness among all of the baryons. From the 
constraint of Eq. (p5|), it follows that more K + 's than K~'s are produced. In other words, 
the difference between K + and K~ abundances provides an indirect measurement of baryon 
number and hence of /it,. 



States, fa) = 
Aot(/Wfi) I_ 

-Btot^s, A*i) B o 



(25) 
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D. Resonance Decays 



For single-generation decays, the resonance contribution to the Wigner function is given 



by P3J25 



r Q~r)a r r°° 
ns-av^P) = H / —g J d A xp L drpTp exp(-r i ar )S ) 





x6 4 (x- 







T /3 

X/3 H P/3 



$/9-a(p^,p)Sj r (^,P i 9) 



(26) 



where the sum in /? is over each decay channel of each resonance that will produce a particle 
a. The quantity §p^ a (pp,p) is the probability density that a resonance with momentum pp 
will decay into a particle a with momentum p. For a two-body decay (3 — > a + X that is 
isotropic in the rest frame of the resonance, 



(Pf3,P) 



-5 E (mp;m a ,m x ) 



PfrP 



where bp is the branching ratio of the particular decay channel and 



(27) 



p (mp;m a ,m x ) = ^—\l[mp 2 ~ {m a + m x ) 2 } [nip 2 - (m a - m x ) 2 } 

ZTTLr 



(28) 



A three-body decay (3 — > a + X + Y can be treated as a two-body decay into a and a 
system with the combined invariant mass M of particles X and K. Since M can vary from 
mx + niy to m/3 — 77i a , it must be integrated over with an appropriate probability density. 



In p4[j26l , this has been shown to be 

p (mp; m a , M)p (M; m x ,m Y ) 



P(M) 



Imx+Zy dMpo(mp; m a , M)p (M; m x ,m Y ) 
Using this probability density, we have for three-body decays 



&P^a+X+Y{Pf3,P) — / dMP(M)<&p^ a+ M 
J mx+my 



(29) 



(30) 



The resonances included which contribute to the charged pion distributions are [27 
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T) — > 7T + + 7T~ + 7T° 

77 — > 7r + + 7r~ + 7 

P — > 7T + 7T 

C<J — > 7T + + 7T~ + 7T° 

C<J — > 7T + + 7T~ 

IT -> K + 7T 

A -> iV + 7T 

£(1385) -> A + 7T 
£(1385) -> S + 7T 

A(1405) -> £ + 7T . (31) 

The X* and A resonances also contribute to the kaon and nucleon distributions, respectively. 

Since the model uses separate isospin, baryon, and strangeness chemical potentials, each 
species of each resonance must be treated separately For example, there are four different 
channels by which negative pions can be produced from delta decay: 

A~ (I = -3/2, B = 1) -> n + n- , 6 = 0.994 

A°(7 = -1/2, B = 1) ->p + 7r - , 6 = 0.994/3 

A~ (/ = -3/2,5 = -1) ->p + 7r _ , 6 = 0.994 

A _ (J = -1/2,5 = -1) ^n + ir- , 6 = 0.994/3. (32) 



III. DETERMINATION OF THE PARAMETERS BY X 2 MINIMIZATION 

A. Construction of \ 2 

A single point i of an experimentally measured one-particle distribution for a particle 
of type a can be characterized by the set (pi, P % a , of at ), where of at is the statistical error 
of the measurement. There are also systematic errors associated with these measurements, 
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which are usually expressed in terms of a percentage of P l a for each particle a. We denote 
the percent systematic error by f a , so that the total error for point i is given by 



<• = \M tat + (f« P «) 2 • (33) 
For each particle a, we construct its contribution Xa 2 to the total x 2 i n the following way: 

Xa (0) - 2^ —2 > (34) 

i u a,i 

where P Q (p, 6) is defined by Eq. (|), 6 is used to represent all of the model parameters (T, 
v t , etc.), and the sum is over all of the measured data points. 

A similar contribution to x 2 can be constructed by comparing two-particle correlation 
data to the corresponding model calculations. The main difference with correlation mea- 
surements is that most of the systematic errors are removed when dividing the two-particle 
distribution measurement in the numerator by the product of the one-particle distribution 
measurements in the denominator. All of the individual one- and two-particle contributions 
to x 2 can be combined into an overall x 2 °f the form 

X [V) - 2^2^ st at2 , (f piV + 2^2^ ~2 ' ^ 

where we explicitly show the dependence of x 2 on the model parameters 9. By varying these 
parameters to minimize x 2 ■, we can find the best fit of the model to the data. Our programs 
can minimize x 2 by using either the Simplex method or the Levenberg-Marquardt method 



B. Confidence in the Overall Fit 

By the central limit theorem, it is well known that the probability distribution of the 
sum of a very large number of small random deviations almost always converges to a normal 
distribution. In the remainder of this paper, we will assume that a sufficiently large number 
of measurements have been taken so that the measurement errors are normally distributed 
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around their true values. With this assumption as well as the assumption of a perfect 
model, the probability density ITj, of obtaining a minimum of x 2 equal to Xmin 2 is given by 
the chi-square function [p5| , p9[] 

n^Xmin 2 ) = 2Y(v/2) ex P(~^ min2 / 2 ) ( ~y~ 1 ' ( 36 ) 

where u is the number of degrees of freedom of the fit. For a model with M adjustable 
parameters fit to N data points, v is simply N — M. Since the above distribution is single- 
peaked with a mean of v and a variance of y/2u, it follows that for a perfect model the 
most probable values of Xmm 2 per degree of freedom are close to one. A value of Xmin 2 /^ 
much larger than one corresponds to a small probability density and consequently leads one 
to question whether the model used is a good one or whether the error bars on the data 
have been underestimated. Conversely, a value of Xmin 2 /^ much less than one seems too 
good to be true and leads one to question whether the error bars on the data have been 
overestimated. 

A more quantitative way to determine the "goodness" of the fit is to integrate U u (x 2 ) 
over x 2 from the Xmin 2 actually found in the fit to infinity. The resulting function V u (Xmin 2 ) 
is the probability that random measurement errors and a perfect model would lead to a 
minimum of x 2 at least as big as the one actually found, Xmm 2 - A fit to a "good" model 
results in a small Xmm 2 and a V u (xmm 2 ) which is greater than some acceptable probability 
(e.g., 5%). A fit to a "poor" model, on the other hand, results in a large Xmin 2 and a 
small value of V u (Xmm 2 )- If, for example, one obtained V u (xmin 2 ) = 10~ 10 , it would be very 
difficult to believe that the large value of Xmin 2 giving rise to that small probability had 
come about purely by way of random measurement errors; it is far more likely that there is 
something seriously wrong with the model. Of course, the choice of which probability should 
actually be used to draw the line between "good" and "bad" models is a purely subjective 
judgement. For our criterion, we choose that probability to be a few percent. 
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C. Estimated Errors in the Parameters 



Assuming the model is determined to be a "good" one, the parameters found at Xmin 2 m &y 
reveal some interesting physics. If, for example, an unexpectedly low freezeout temperature 
is discovered in a fit, it is necessary to know how confident one can be in that particular low 
value. In other words, we need to estimate the possible error in that fitted parameter due to 
random measurement errors in the experiment. More generally, one would like to determine 
the region in M-dimensional parameter space which has a high probability of containing the 
underlying true parameter values. For example, we would like to be able to say "there is a 
99% chance that the true parameter values fall within such and such region of parameter 
space." 

For normally distributed measurement errors, the relevant region is one consisting of all 
combinations of parameters 9 which would lead to a value of x 2 ($) l ess than x min 2 + A, where 
A is some constant number. The confidence level "CL" for a specific A and a specific number 
of parameters M is given by the integral of an M-degree-of- freedom chi-square distribution 

US 

A 





cl = J Q A d x 2 n M (x 2 ) (37) 



For example, if we choose a region defined by A = 21.666 for a model with M = 9, Eq. ( |3"TD 
tells us that CL = 0.99. This means that we have 99% confidence that the true parameter 
set lies among all possible sets which result in a x 2 (#) less than Xmm 2 + 21.666. By inverting 
Eq. (|3~7"D , it is possible to find the appropriate A for M parameters which leads to any desired 
confidence level. 

Once A has been determined in this way, the desired region in parameter space is found 
by making a Taylor expansion of x 2 {@) about its minimum at 9: 



1 M - - 8 2 y 2 

X 2 (0) = Xmin 2 + ^ E ~ a )(9 b - 9J v 



1 a r=x 'd9 a d9 h 



(38) 



Notice that the gradient terms disappear because we are expanding around a minimum. If 
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we assume that higher-order terms are also negligible close to the minimum, then the desired 
region of parameter space is just the M-dimensional hyperellipsoidal volume defined by 



M 



E % - a )D. 



ab\Vb 



Oh) < A 



(39) 



a,b=l 



where the curvature matrix D ab is just one-half of the second derivative matrix of x 2 . 

From Eq. (|35|), it is apparent that the exact curvature matrix D ab involves terms of the 
form 



and 



af at 2 + {f a piy 



Oi 



d 2 P a ( Pl ,9) 

de a de b 



d 2 Cp(q i} K i} < 



d9 a d9 b 

For a good model, P % a — P a and C % a — C a are nothing more than the random measurement 
errors at the point i. Since these errors can have either sign and should in general be 
uncorrelated with the model, terms of the above form tend to cancel out when summed over 
i. By dropping these terms, we obtain the approximate curvature matrix that is actually 
used in the model: 

i fdP a ( Pi ,6)\ (dP a (pi,ey 



a i 

+ EE 



af at2 + (f a p^y 

1 (dCpiq^Ki 



P i 



Of 



d9 a J \ d9 b 

\ de h ) 



(40) 



The best and most complete way to specify the parameter confidence region is to calculate 
the above curvature matrix and then employ Eq. fl3~9|) . However, in order to get a quick idea 
of the size of the region, it is useful to provide one- dimensional error estimates on each 
parameter. One way to do this is by determining the largest and smallest values of each 
parameter that can be obtained on the hyperellipsoid. For example, suppose we had only 
two parameters, T and v t , and we determined that the confidence region was the interior of 
the ellipse shown in Fig. 3. The one-dimensional error estimates 5T and Svt shown in the 
figure would give one a rough idea of the size of the ellipse. One must be careful in using 
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these one-dimensional simplifications, however, since for example the point (T + ST,v t + 5v t ) 
lies outside the confidence region. Mathematical determination of the one-dimensional error 
estimates is derived in Appendix A for the adjusted parameters and in Appendix B for 
additional calculated quantities of physical interest. 



IV. FITTING TO E-802 DATA 



A. Description of the Data 



Having developed our model, we used it to fit meson data from central Si + Au collisions 
at piab/A = 14.6 GeV/c measured in experiment E-802 at the AGS As mentioned 

previously, we did not attempt to fit proton or deuteron data due to the difficulty of disen- 
tangling spectators from participants. The invariant one-particle multiplicity distributions 
for 7r + , 7r~, K + , and K~ were obtained from an on-line database |7[] and then organized into 
input files in which each line contained the information 



y P ,i pm pl a 



i ^stat 

i 



These data were presented in graphical form in ||, where it was estimated that systematic 
errors for all of the particles except A~'s were 10-15%, while those for the A~'s may have 
been as much as 20%. For our fits, we used constant systematic errors of 15% for pions and 
A+'s and 20% for A"'s. 

In addition to these data, we obtained preliminary three-dimensional correlation data 
for both 7r + and K + ||. For the correlation data, the momenta of the identical particles 
were measured relative to a fixed frame at yi a b = 1.25. Each two-particle event was then 
three-dimensionally binned according to the components q z , q out , and g si d e of its momen- 
tum difference. We use here the standard notation of z denoting the beam direction, "out" 
denoting the direction parallel to the component of the average momentum K which is per- 
pendicular to the beam, and "side" denoting the remaining transverse direction [3C|. The 



convention was used that g out = Pi,out — £>2,out is always positive [pl] . Since this convention 



identifies which particle is "particle 1" and which is "particle 2," the remaining two momen- 
tum differences, q z = p± tZ — p 2 , z and g si d e = Pi, s ide — £>2,sidc, can and did take either sign. In 
other words, there were bins with negative values of q z and q S id c as well as those with positive 
values. For each bin in q, the average values of Y — \{yi + y-i) and K t = |(pi iOU t + P2,out) 
for the pairs in that bin were calculated and recorded. Thus, in the correlation input files, 
each line contained the information 

00 » \Kt)i Qz,i Qout,i <lside,i C % Oi , 

where ()$ represents the average value over the bin i. 

Due to the very large number of correlation data points, not all of the data were used. 
In particular, for n + we used only points for which the at the center of the bin satisfied 
the inequality 

f^v + (?^y + f^v < i > ( 4i ) 



V2oo/ uooy v 100/ 

with the components of q measured in MeV/c. Similarly, for K + we used only data points 
satisfying 

Our rationale for omitting some of the data points was that most of the physically important 
information is contained at low momentum differences. For momentum differences greater 
than those specified above, the model predicted a correlation function very close to one, 
while the actual data fluctuated wildly about this value with huge error bars. Furthermore, 
we found that the use of different momentum-difference cutoffs did not significantly affect 
our results. 

Combining both the one-particle distributions and the two-particle correlations, we used 
a total of 1416 data points in our fits. Since there are nine adjustable parameters, there are 
1407 degrees of freedom. 
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B. Results of the Fit 



The best-fit parameter values and their estimated errors (see Appendix A) at a 99% 
confidence level for the minimum \ 2 found are listed in Table I. Since the 99%-confidence 
hyperellipsoid overlaps an unphysical region in which a t < — 1, the lower error estimate on at 
has been reduced to reflect only the physical region. The curvature matrix D and its inverse, 
the covariance matrix D -1 , are presented in Tables II and III, respectively. Having found 
the minimum, it is possible to calculate a number of other related quantities of physical 
interest. Some of these are constrained parameters such as /i s , while others are calculated 
quantities such as the maximum longitudinal expansion velocity V£ = tanh(^o) achieved by 
the source in its center-of-mass frame. One of the most interesting of these quantities is the 
local baryon density at freezeout. Appendix B describes how each of the related quantities 
is calculated, while Table IV shows all of these quantities with their estimated errors. 

The "goodness" of the fit can be seen by looking at Fig. 4, which shows a plot of the 
probability density U u of Eq. (|36"D for 1407 degrees of freedom versus Xmin 2 - The high 
probability density of the resulting Xmm 2 = 1484.6 is good evidence that there is nothing 
seriously wrong with the model. In fact, integration of the shaded region determines that 
there is a 7.4% chance of obtaining a value of Xmm 2 at least that large when fitting an 
absolutely perfect model to the data. To get an idea of how much each data set contributed 
to the overall Xmin 2 , we have broken down the individual contributions in Table V. Obviously, 
all of the data sets are fit reasonably well, although the \ 2 P er degree of freedom for the 
negative-pion one-particle distribution is slightly higher than the rest. 

The goodness of fit for the one-particle distributions can also be seen qualitatively by 
looking at direct comparisons of the model to the data. Figures 5-8 show the theoretical and 
experimental meson invariant one-particle multiplicity distributions as functions of m t — m 
for various rapidities. As expected from the low value of Xmm 2 , the overall agreement looks 
excellent. 

Although proton data were not used for the main fit, once the best parameters have been 
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found, it is possible to calculate the proton distribution and compare it with the experimental 
data. Figure 9 shows that the model predictions agree with the proton data moderately well 
for rapidities of 1.3 or greater. For lower rapidities, however, the data show far more low-p t 
protons than are predicted by the model. In our picture, we consider these protons to be 
target spectators which may have interacted somewhat, but not enough to be considered 
part of the hydrodynamical system. It should also be noted that since our model does not 
distinguish between a deuteron and a separate proton and neutron, some of the excess when 
the model overpredicts the proton data may be due to leaving out deuteron coalescence. 

As mentioned before, the main problem with including proton data in the fits is figuring 
out how to eliminate contamination by spectators. Nevertheless, in order to get some idea 
of how the inclusion of protons might affect our results, we made four different fits in which 
we included all proton data in the following rapidity ranges: 1.5 < y p < 1.9, 1.3 < y p < 1.9, 
1-5 < y p < 2.1, and 1.3 < y p < 2.1. These fits all give extremely similar results, so we 
will discuss only the last case. The 190 proton data points for this case bring the number 
of degrees of freedom up to 1597. The minimum \ 2 of 1730.4 (x 2 /V = 1-084) is obtained 
with the following parameters: T = 95.8 ± 3.9 MeV, fi h /T = 5.30 ± 0.28, X n = 0.66 ± 0.11, 
R = 7.8±1.5 fm, v t = 0.640 ±0.034 c, a t = -0.87 t° f 3 , y s = 1.320 ±0.072, r] = 1.48 ±0.14, 
and Tf = 8.1 ± 2.1 fm/c. The central values themselves for all of the parameters except 
fib/T lie within the individual 99%-confidence intervals of the original fit (see Table I). Even 
for fib/T, the 99%-confidence interval for the fit with protons overlaps that for the original 
fit. In other words, inclusion of the proton data changes the best-fit parameter values only 
within their stated uncertainties. It does, however, lead to a significant reduction in the 
calculated number of projectile participants. For the fit with protons included, the number 
of projectile participants is reduced to 17.7 ± 1.4, which is to be compared to 26.1 l®;® found 
in Table IV and to 28 nucleons in a 28 Si nucleus. If one were to take the proton fit seriously, 
one might wonder how a central Si ± Au collision could possibly give rise to 10 projectile 
spectators. It is our view that the unresolved issues of coalescence and spectator separation 
in our model make it better to simply neglect the protons altogether, just as we did in our 
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original fit. For the remainder of the paper, we will always refer only to that original fit. 

Since the preliminary correlation data that was used has not yet been published, we show 
here in Figs. 10 and 11 (g 2 ,g ut) projections of the correlation functions calculated by the 
model. Notice that whereas the kaon correlation function intercepts the q = axis at 2, the 
pion correlation function intercepts the axis at 1.65, corresponding to a value of = 0.65. 
Notice also that since nonvanishing values of both K t and Y — y s are used for the plots, 
the effects of the "out-long" cross term |p2"l,[32"l, j33t can be seen (especially in the kaons) as a 



slight twisting in the major axes of the correlation function. 

Data from central Si + Au collisions at the AGS have also been compared to a thermal 



model in [34J. There it was argued that a freezeout temperature of 120-140 MeV was 
consistent with these data. We, on the other hand, have found a much lower temperature 
of 92.9 ± 4.4 MeV (see Table I). To explain this significant discrepancy, we would like to 



point out a few differences between our approach and that of |34| . First of all, they assumed 



that transverse and longitudinal flow are completely separable. This led them to compare 
a thermal model that had been integrated over all rapidities [^TJ with data from a single 
mid-rapidity bin. In contrast, we make no such assumption. Secondly, the model used in 
34J never specifies the size or shape of the freezeout hypersurface. In addition to preventing 



comparisons to two-particle correlations, this ambiguity forces them to multiply each of 
their one-particle distributions by an arbitrary normalization factor before comparing to the 
normalized data of 0. Not only does our unambiguous parametrization of the freezeout 
hypersurface allow us to compare to two-particle correlations, it also allows us to see the 
significant effects that different temperatures have on the absolute normalizations of one- 
particle distributions. Another significant difference between the two approaches is that 
in |34| only two points in temperature were studied, whereas in our approach the whole 
nine- dimensional parameter space is explored, resulting in the absolute minimum of x 2 . 

In order to see how much worse higher temperature fits would be, we made a number of 
runs at various fixed temperatures, allowing all of the other parameters to vary. The results 
are given in Table VI and in Fig. 12, a plot of the minimum \ 2 & t a fixed temperature vs. 
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that temperature. As Table VI shows, the value of x 2 f° r T = 120 MeV is 2025.4. Again by 
integrating Tl u of Eq. ([36]), we can determine that the probability of a perfect model resulting 
in a x 2 a t least as large as 2025.4 is the incredibly small value 5.1 x 10~ 25 . Above T = 129 
MeV, the minimum x 2 solution switches to a different branch. This high-temperature branch 
is actually unphysical, as can be seen by examining Table VI and noting that it is impossible 
for the system to expand to the large transverse radius R in the infinitesimally small time 

tl = Tfy/1 + a t . 

Another result of the model is the size and shape of the freezeout hypersurface. Since 
at for the best fit is negative, freezeout begins at time t — t\ — 3.1 fm/c at z = and 
p = R = 8.0 fm, and takes 5.1 fm/c to reach the center of the source at z = p = at 
time t = ti = Tf = 8.2 fm/c. Freezeout along the symmetry axis then occurs at a constant 
proper time, finally ending at source-frame time = Tf cosh 770 = 18.8 fm/c. As mentioned 
previously, Figs. 1 and 2 pictorially show the freezeout process for these parameter values. 

V. FUTURE ISSUES 

The actual "freezeout" process taking place in these collisions is undoubtedly far more 
complicated than in our model. Azimuthal symmetry may be broken, the local temperature 
and chemical potentials may have some spacetime dependence, the expansion flow veloc- 
ity may be neither boost-invariant nor linear in p, the hypersurface may have a different 
shape and/or some four- dimensional fuzziness, kaons may freeze out before pions, chemical 
equilibrium may not be fully achieved, etc. One may even question whether equilibrium hy- 
drodynamical concepts are valid at all. Although this paper does not definitively settle these 
questions, the remarkable agreement between theory and experiment suggests that our real- 
istic nine-parameter expanding source model nevertheless provides a very good description 
of the most important physics taking place at freezeout. 

One parameter which definitely needs to be better understood is the incoherence pa- 
rameter A^. Does the fact that it is less than one mean that a significant number of pions 
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are being produced coherently, or could the reduced intercept instead be largely an artifact 
arising from the way the correlation function was determined experimentally ||35|| ? 

We hope that our model will be used in the future to systematically analyze the de- 
pendence of the freezeout quantities upon bombarding energy and the sizes of the colliding 
nuclei. A sharp discontinuity in one or more of these quantities could be a signal of quark- 
gluon-plasma formation. 
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APPENDIX A: ONE-DIMENSIONAL ERROR PROJECTIONS 

The simplest way to determine the largest and smallest values attained by parameter 9 a 
on the hyperellipsoid defined by Eq. (^) is to use a Lagrange multiplier £ a . We begin by 
finding the maximum (or minimum) of the quantity 

M 

6 a — £,a ^2 (@b — 9b)D bc (9 c — 6 C ) . 

b,c=l 

By differentiating with respect to 9d, we find the coordinates 8 e of the extrema as a function 

«.-».=^(b-'L. < A1 > 

where the subscripts a are not summed over. To impose the constraint that the solution 
lies on the hyperellipsoid, we must pick a £ a such that the equality of Eq. is satisfied. 
Plugging Eq. ( |Al| ) into Eq. (|3^) and solving for £ a , we find 
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A 



aa 



(A2) 



where again the indices a are not summed over. Inserting Eq. ( |A2|) into Eq. (|Al|), we find the 



values of all of the parameters 6 e corresponding to each extremum of 6 a on the hyperellipsoid. 
In particular, the one- dimensional error estimate on 9 a is just proportional to the square 
root of the aath element of the covariance matrix (the inverse of the curvature matrix): 



By inverting the curvature matrix to get the covariance matrix, one can then just read off the 
one-dimensional error estimate on each parameter by taking the square root of the product 
of the appropriate diagonal element times A. 

APPENDIX B: ADDITIONAL CALCULATED QUANTITIES 

Here we list some additional physical quantities of interest which can be calculated from 
the nine parameters. From Eqs. (|7|) and (0), it can be seen that the maximum longitudinal 
velocity achieved by the source is given by v# = tanh^o- Also, in Sec. II. B the times that 
freezeout begins, reaches the center of the source, and ends were shown to be ti = + a t , 
t 2 = Tf, and t 3 = n cosh?7o, respectively. The maximum longitudinal extension of the source 
is given by z 3 = Tf sinh i] , while the duration of freezeout at z = is given by At = t 2 — 1\. 
The total baryon number B tot of the source is given by the denominator of the second 
equation in (pip. Section II. C also explains how the chemical potentials fi s and /ii are found. 
The numbers of projectile and target participants can be deduced from Eq. ( p3|) and are 




(A3) 



given by 



5, 



B tot sinh(y s - y tar ) 



proj — 



sinh(y s - y tar ) + smh(y proi - y B ) 
B tot smh(y pioi - y s ) 



B, 



'tar — 



sinh(y proj - y s ) + sinh(y s 

— 2/tarJ 



(Bl) 



The local density of particles of type a is given by the integral 
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n a {p) = —. — / ^ / dtS(x,p) 



7 P cosh 7] J E 



tP 



IP k = l 

B?Jl + a t {p/RY °{ T ) l { RT 



[B2) 



where Ki and Jj are modified Bessel functions of order i. Due to the boost invariance 
assumed in everything but the spatial limits of the model, n a is a function of only p and not 
of r). The local baryon density is just given by 

n M = ^2B a n a (p) . (B3) 

a 

To calculate error estimates on these quantities, we could in principle use a more general 
form of the Lagrange-multiplier method introduced in Appendix A. We have found, however, 
that a quicker and more reliable method is to first find M — 1 new variables which can 
parametrize just the surface of the hyperellipsoid, express the quantities as functions of 
these new variables, and then find the extrema of these functions. We begin this process by 
numerically finding the unitary matrix U which transforms D into the diagonal matrix D, 
namely 

D = U- l DU. (B4) 
Next we use the (diagonal) elements of D to define the vector 



/ f) M 



where there is no summation over the index a. Using these variables, we can see that the 
equality of Eq. (p9|) reduces to the equation of the surface of a sphere in M dimensions, 
namely 



E^a 2 = L (B6) 



a=l 



Since the surface is (M — l)-dimensional, M — 1 angles are sufficient to identify any point 
on it. We now redefine the ip a in terms of the angles (f) a through 
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^2 



= COs(0i) 

= sin(0i) cos(0 2 ) 



ipM-i = sin(0i) sin(0 2 ) • • •sin(0 M _ 2 ) cos(0 M -i) 
V>m = sin(0i) sin(0 2 ) • • • sin(0 M „ 2 ) sin(0 A/ „i) . (B7) 

Since the inverse of Eq. ( [Bop tells us 



nr M 

Oa = 0a + \ T -Y. U ^b, (B8) 

we now have the M parameters 9 a expressed in terms of the M — 1 parameters <p a . As 
mentioned previously, the extrema on the hyperellipsoid for any function f(9(4>)) can be 
found simply by allowing the <p a to vary freely. 
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FIGURES 

FIG. 1. Freezeout hypersurface, which specifies the positions in spacetime where the expanding 
hydrodynamical fluid is converted into a collection of noninteracting, free-streaming hadrons. The 
parameters used are R = 8.0 fm, Tf = 8.2 fm/c, ryo = 1-47, and a t = —0.86. 

FIG. 2. Seven snapshots at equal spacings in the source- frame time t of the part of the source 
of Fig. |l] that is freezing out (inner surfaces) or has not yet frozen out (end caps). To draw the 
end caps at each time t, we used v t = 0.683 c and assumed that none of the fluid was accelerated 
prior to reaching the freezeout hypersurface. 

FIG. 3. An example 99%-confidence ellipse for a model with two parameters showing how 
one-dimensional errors on those parameters are estimated. 



FIG. 4. The probability density II„ of Eq. (36) for obtaining a particular minimum of x from 



a fit with 1407 degrees of freedom v to a perfect model. There is a 7.4% chance that a perfect 
model would have given rise to a Xmin 2 at least as large as the one actually found, 1484.6. 

FIG. 5. Comparison between model predictions and experimental data [6,7] for the invariant 
7r + one-particle multiplicity distribution Ed?N/dp 3 = l/(27rmt) d 2 N/dy p dmt as a function of y p 
and mt — m. For visual separation, the results at a particular particle rapidity y p relative to the 
laboratory frame are scaled by the factor 10 5 (°- 7 ~yp) . Although not always distinguishable on the 
scale of the graphs, statistical errors are given by the inner error bars, and total errors are given 
by the outer error bars in Figs. 5-9. 

FIG. 6. Comparison between model and experimental ir~ distributions, scaled as in Fig. ||. 
FIG. 7. Comparison between model and experimental K + distributions, scaled as in Fig. [|. 
FIG. 8. Comparison between model and experimental K~ distributions, scaled as in Fig. ||. 
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FIG. 9. Comparison between model and experimental proton distributions, scaled as in Fig. ||. 
Note that these proton data were not included in the fit. 

FIG. 10. Dependence of the predicted ir + two-particle correlation function C upon the longi- 
tudinal and "out" momentum differences, for fixed values of the other three quantities upon which 
C depends. 

FIG. 11. Dependence of the predicted K + two-particle correlation function C upon the longi- 
tudinal and "out" momentum differences, for fixed values of the other three quantities upon which 
C depends. 

FIG. 12. Minima of x 2 f° r fixed values of T, when all other parameters are allowed to vary. 
Solid circles connected by a solid line identify one branch which contains some physically relevant 
solutions near the minimum, while open circles connected by a dashed line represent a different 
branch containing only unphysical solutions (see Table VI). 
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TABLES 

TABLE I. Nine adjusted source freezeout parameters. 





Value and uncertainty 


Parameter 


at 99% confidence 



Nuclear temperature T 
Baryon chemical potential [ih/T 
Pion incoherence fraction 
Transverse freezeout radius R 
Transverse freezeout velocity v t 
Transverse freezeout coefficient a t a 
Source rapidity y s 
Longitudinal spacetime rapidity rjo 
Longitudinal freezeout proper time Tf 
a Physically, at cannot be less than — 1 , since that value corresponds to freezeout beginning imme- 
diately at t = 0. 



92.9 ± 4.4 MeV 
5.97 ± 0.56 
0.65 ± 0.11 
8.0 ± 1.6 fm 

0.683 ± 0.048 c 

-0.86 ±°f 4 
1.355 ± 0.066 
1.47 ± 0.13 

8.2 ± 2.2 fm/c 
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TABLE II. Curvature matrix D. Derivatives with respect to the parameters are ordered in 
rows and columns in the same way as for the parameters in Table I. The units of each element are 
given by the inverse of the units of the associated row and column parameters (e.g., the units of 
D u are l/(MeV-fm)). 



65.72 


117.4 


-11.16 


175.6 


1859 


-54.69 


606.7 


524.4 


140.5 


117.4 


393.1 


-9.373 


295.8 


1717 


-36.53 


867.3 


740.1 


222.9 


-11.16 


-9.373 


2017 


-127.6 


368.7 


100.0 


-110.7 


-124.2 


-82.88 


175.6 


295.8 


-127.6 


509.4 


4641 


-169.9 


1761 


1501 


395.3 


1859 


1717 


368.7 


4641 


87117 


-2530 


13964 


14133 


4107 


-54.69 


-36.53 


100.0 


-169.9 


-2530 


313.5 


-678.5 


-580.9 


-123.0 


606.7 


867.3 


-110.7 


1761 


13964 


-678.5 


19418 


10296 


1315 


524.4 


740.1 


-124.2 


1501 


14133 


-580.9 


10296 


7643 


1158 


140.5 


222.9 


-82.88 


395.3 


4107 


-123.0 


1315 


1158 


319.0 



TABLE III. Covariance matrix D . The ordering of the rows and columns is the same as in 
Table II. The units of each element are 10 4 times the inverse of the units in Table II (e.g., the 
units of (D-^u/W 4 are MeV-fm). 



8860 


-891.5 


-53.42 


-1180 


-61.21 


-121.3 


-11.63 


-29.23 


-935.8 


-891.5 


143.9 


2.678 


76.55 


7.407 


5.523 


1.195 


4.962 


81.77 


-53.42 


2.678 


5.827 


6.837 


-0.04256 


-1.918 


0.04292 


-0.3064 


15.44 


-1180 


76.55 


6.837 


1244 


24.25 


146.0 


1.241 


-8.943 


-1302 


= -61.21 


7.407 


-0.04256 


24.25 


1.080 


4.091 


0.1796 


0.04903 


-21.53 


-121.3 


5.523 


-1.918 


146.0 


4.091 


61.51 


0.1794 


0.9064 


-164.8 


-11.63 


1.195 


0.04292 


1.241 


0.1796 


0.1794 


1.998 


-3.112 


3.576 


-29.23 


4.962 


-0.3064 


-8.943 


0.04903 


0.9064 


-3.112 


8.359 


2.626 


-935.8 


81.77 


15.44 


-1302 


-21.53 


-164.8 


3.576 


2.626 


2193 
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TABLE IV. Additional calculated physical quantities. 





Value and uncertainty 


Quantity 


at 99% confidence 



Source velocity v s 
Longitudinal velocity vi 
Longitudinal freezeout radius z% 
Beginning freezeout time t\ 
Freezeout time ti at source center 
Final freezeout time ts 
Freezeout width At in proper time a 
Baryon chemical potential /j,^ 
Strangeness chemical potential /j, s 
Isospin chemical potential //; 

Number -B pro j of baryons originating from projectile 
Number B tar of baryons originating from target 
Total number B to t of baryons in source 
Baryon density n\ at beginning of freezeout b 
Baryon density n s along symmetry axis 
a Calculated under the additional assumption that the exterior matter at z = that freezes out 
first has been moving with constant transverse velocity vt from time t = until time t±. 
b The upper limit of oo for this quantity arises because the beginning freezeout time t\ could be 
zero, at which time the shape is an infmitesimaHy thin disk. 
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U.O/D _ .016 C 

0.900 t° fi c 
16.9 tH fm 
3.1 ±U We 
8.2 ± 2.2 fm/c 
18.8 fm/c 
5.9 tti fm/c 
554 tit MeV 
75 t\l MeV 
-5.3 tl:? MeV 
26.1 l|;f 

c; 7 +20 
°' -15 

83 tjf 
0.057 t - 032 fm- 3 

0.0222 ±8:88^ fm- 3 



TABLE V. Individual contributions to x 2 ■ 



Type of data iY data v & X 2 X 2 /v 



7r + one-particle 


231 


229.5 


238.0 


1.037 


ir~ one-particle 


239 


237.5 


266.2 


1.121 


K + one-particle 


137 


136.1 


140.6 


1.033 


K~ one-particle 


49 


48.7 


51.2 


1.051 


7r + correlation 


464 


461.1 


498.0 


1.080 


K + correlation 


296 


294.1 


290.7 


0.988 


Total 


1416 


1407.0 


1484.6 


1.055 



a In determining the number of degrees of freedom for the individual contributions, we have allocated 
the nine adjustable parameters among each type of data in proportion to the number of data points. 
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TABLE VI. Eight adjusted parameters of best fits at fixed values of temperature T plus cal- 
culated -Bproj- Physically relevant solutions correspond only to a limited region surrounding the 
absolute minimum at T = 92.9 MeV. For values of T below 82 MeV, the calculated lower limit 
on -Bproj exceeds the number of nucleons in the projectile. Solutions below the horizontal line 



corresp 


ond to a new 


r branch 


which is 


unphy 


sical for the reason 


mentioned 


in the 


text. 




x 2 


T 






R 


v t 




y s 


% 


Tf 


-Bproj 




(MeV) 






(fm) 


(c) 








Cfm/cl 




3280.7 


50 


16.4 


1.02 


21.1 


0.971 


0.232 


1.467 


1.75 


10.2 


566 


2557.2 


60 


12.7 


0.868 


14.2 


D Q9D 




1.428 


1.70 


10.8 


204 


2043.2 


70 


9.98 


0.777 


11.1 




— n R7i 

u.u / ± 


1.394 


1.64 


10.5 


88.9 


1678.9 


80 


7.76 


0.723 


9.64 


PI 778 


— n 7^s 


1.375 


1.55 


9.65 


41.9 


1494.6 


90 


6.28 


0.667 


8.39 


n 70^ 


— n 8i 8 


1.361 


1.49 


8.47 


27.8 


1534.4 


100 


5.35 


0.613 


7.07 


vj.vj^j ± 


— QQQ 


1.341 


1.46 


7.55 


23.6 


1739.1 


110 


4.68 


0.544 


6.13 


n ^fi8 

u.uuo 


— D QQQ 


1.315 


1.47 


6.44 


21.6 


2025.4 


120 


4.17 


0.486 


5.24 


VJ .UVJ£t 


— QQQ 


1.284 


1.49 


5.69 


20.1 


2323.9 


130 


3.81 


0.473 


3.42 


vj.ooo 


— QQQ 


1.266 


1.54 


7.40 


19.2 


2437.8 


135 


3.70 


0.545 


2.37 


VJ.Zj<D± 


— QQQ 


1.268 


1.57 


10.3 


19.3 


2304.5 


125 


7.17 


0.530 


13.5 


879 


QQQ 


1.380 


1.47 


0.133 


72.5 


2294.2 


130 


7.05 


0.530 


13.3 


u.ouu 


— QQQ 


1.379 


1.47 


0.102 


68.2 


2273.1 


140 


6.85 


0.529 


13.0 


U.004 


— QQQ 


1.375 


1.45 


0.0630 


60.3 


2252.2 


150 


6.69 


0.528 


12.7 


U.OOj 


— QQQ 

VJ. 


1.371 


1.43 


0.0407 


53.2 


2230.9 


160 


6.60 


0.528 


12.4 


0.823 


-0.999 


1.268 


1.41 


0.0269 


45.8 


2209.3 


170 


6.54 


0.530 


12.1 


0.804 


-0.999 


1.365 


1.39 


0.0183 


37.5 


2190.9 


180 


6.41 


0.532 


11.9 


0.783 


-0.999 


1.364 


1.36 


0.0139 


31.1 


2182.1 


190 


6.20 


0.533 


11.6 


0.760 


-0.999 


1.361 


1.35 


0.0116 


28.4 


2185.2 


200 


5.96 


0.536 


11.4 


0.738 


-0.999 


1.359 


1.33 


0.0103 


27.7 


2199.9 


210 


5.76 


0.537 


11.1 


0.713 


-0.999 


1.353 


1.32 


0.00905 


26.4 


2228.1 


220 


5.57 


0.537 


10.9 


0.690 


-0.999 


1.347 


1.32 


0.00805 


27.0 


2271.9 


230 


5.28 


0.542 


10.7 


0.666 


-0.999 


1.333 


1.33 


0.00765 


26.7 


2335.3 


240 


4.99 


0.544 


10.5 


0.637 


-0.999 


1.325 


1.34 


0.00756 


28.1 


2407.0 


250 


4.73 


0.545 


10.2 


0.605 


-0.999 


1.319 


1.34 


0.00748 


28.7 
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Si + Au, central 7% 
P\JA= 14.6GeV/c 




20 



10 t 

(fm/c) 



20 



10 * 



Source freezeout 
Si + Au 5 central 7% 
P\*JA= 14.6GeV/c 




Time t (fm/c) = 
3.10 



5.65 



8.20 



10.75 



$ . 

a 13.30 

I 15.85 
I 18.40 

Figure 2 



39 




1 > 

T 



Figure 3 



0.008 




1200 1300 1400 1500 1600 

Minimum of % 2 



Figure 4 




0.0 0.5 1.0 1.5 2.0 
Transverse Kinetic Energy (GeV) 



Figure 5 
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Figure 6 
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Figure 7 
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Figure 8 
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Figure 9 
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Figure 12 



